/*
NOTES TO DATA ANALYSIS:
	1. Treatment effect regressions (each regression should control for strata variables)
		- by pooled treatment
		- by 5 frames
		- by 2 timings
		- by 10 arms
*/

clear all
set more off
capture log close

********************************************************************************

log using "$logfiles/analysis_key_itt_`c(current_date)'.log", replace

use "$covidclean/smscovid_clean.dta", clear

***Keep relevant observations: consent taken and above 18 years of age, non-missing treatment arms
	keep if consent==1 & age>=18
	keep if ~missing(treatment_arm)
	count //3964

***HOUSEKEEPING - BEHAVIOR TREATMENTS, OUTCOMES, ROW TITLES

	***Create dummies with treatment arms interacted with behavior
		foreach v in $treat_pool $treat_frames $treat_timings $treat_arms {
			gen `v'_sd = `v'*behavior_sd
			gen `v'_hw = `v'*behavior_hw
		}
		
	***Outcomes
		global keysdhw know_sd act_sd know_hw act_hw
	
	***Row titles
		global rtitle_pooled "Treatment - SD \ "" \ "" \ Treatment - HW \ "" \ "" "
		global rtitle_frames "SD - Neutral \ "" \ "" \ SD - Public Gain \ "" \ "" \ SD - Public Loss \ "" \ "" \ SD - Private Gain \ "" \ "" \ SD - Private Loss \ "" \ "" \ HW - Neutral \ "" \ "" \ HW - Public Gain \ "" \ "" \ HW - Public Loss \ "" \ "" \ HW - Private Gain \ "" \ "" \ HW - Private Loss"
		global rtitle_timings "SD - 2xmorn. \ "" \ "" \ SD - Morn./Eve. \ "" \ "" \ HW - 2xmorn. \ "" \ "" \ HW - Morn./Eve."
		global rtitle_arms "SD - Neutral (2xmorn.) \ "" \ "" \ SD - Public gain (2xmorn.) \ "" \ "" \ SD - Public loss (2xmorn.) \ "" \ "" \ SD - Private gain (2xmorn.) \ "" \ "" \ SD - Private loss (2xmorn.) \ "" \ "" \ SD - Neutral (morn./eve.) \ "" \ "" \ SD - Public gain (morn./eve.) \ "" \ "" \ SD - Public loss (morn./eve.) \ "" \ "" \ SD - Private gain (morn./eve.) \ "" \ "" \ SD - Private loss (morn./eve.) \ "" \ "" \ HW - Neutral (2xmorn.) \ "" \ "" \ HW - Public gain (2xmorn.) \ "" \  "" \ HW - Public loss (2xmorn.) \ "" \ "" \ HW - Private gain (2xmorn.) \ "" \ "" \ HW - Private loss (2xmorn.) \ "" \ "" \ HW - Neutral (morn./eve.) \ "" \ "" \ HW - Public gain (morn./eve.) \ "" \ "" \ HW - Public loss (morn./eve.) \ "" \ "" \ HW - Private gain (morn./eve.) \ "" \ "" \ HW - Private loss (morn./eve.)"

********************************************************************************
***MAIN ITT REGRESSIONS
********************************************************************************
	
		qui foreach v of varlist $keysdhw {
			
			********************************************************************
			***BY POOLED TREATMENT
			********************************************************************
			
				*Control mean
				mean `v' if treatment_pooled==0 & $sample
				mat mean = r(table)
				local c_mean : di %9.2f mean[1,1]

				*Regression
				noi reghdfe `v' $treat_pool_sd $treat_pool_hw if $sample, a($studycontrols $covariates) vce(robust)

				mat b = r(table)
				local tcount: word count $treat_pool_sd $treat_pool_hw				
				forval i = 1/`tcount' {
					local b_`i' = b[1,`i'] //beta
					local se_`i' = b[2,`i'] //asymptotic s.e.
					local a_p_`i' = b[4,`i'] //asymptotic p-value
					local p_`i' = 0 //exact p-value
				}
				
				local n : di %9.0fc e(N)
				local r2_a : di %9.2f e(r2_a)
				
				*Randomization inference using RI dataset created
				preserve
    				forval batch = 1/$nbatch {
					   use if batch == `batch' using "$rianalysis/ri_data_final", clear
    				   noi di "Running ITT-pooled for `v' for batch - `batch' at time: `c(current_time)'"
    				   reghdfe `v' $treat_pool_sd $treat_pool_hw if $sample, a($studycontrols $covariates) vce(robust) 
    				   mat ri_b = r(table)
    				   forval i = 1/`tcount' {
            			   local ri_b_`i' = ri_b[1,`i']
            			   *One-tailed comparison
				           if `ri_b_`i'' >= `b_`i'' local p_`i' = `p_`i'' + 1 
					    }  
				    } 
				restore
				
				*Compute Exact p-values 
				forval i = 1/`tcount' {
					local p_`i' = `p_`i''/$nbatch
				}
				
				*Imputing coefficients, asymptotic and exact p-values in the output matrix
				mat `v'_p = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat `v'_p[`i',1] = round(`b_`i'',0.001)
					mat `v'_p[`i',2] = round(`se_`i'',0.001)
					mat `v'_p[`i',3] = round(`p_`i'',0.001)
				}
								
				*Statistical significance
				mat stars = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat stars[`i',2] = 0
					mat stars[`i',3] = 0
					if (`a_p_`i'' < 0.1)  mat stars[`i',1] = 1
					else mat stars[`i',1] = 0
					if (`a_p_`i'' < 0.05) mat stars[`i',1] = 2 
					if (`a_p_`i'' < 0.01) mat stars[`i',1] = 3 
				}
				
				*Save matrix
				matsave `v'_p, replace saving path($tables)
				*matload `v'_p, path($tables) saving
				
				*Output 
				frmttable, statmat(`v'_p) store(`v'_p) substat(2) annotate(stars) asymbol(*,**,***) rtitle($rtitle_pooled) /// 
				addrows("Adjusted $ R^2$", `r2_a' \ N, "`n'" \ "Control Mean", `c_mean') sdec(3,3,3)
						
			********************************************************************
			***BY 5 FRAMES
			********************************************************************
				
				*Control mean
				mean `v' if treatment_pooled==0 & $sample
				mat mean = r(table)
				local c_mean : di %9.2f mean[1,1]
				
				*Regression 
				noi reghdfe `v' $treat_frames_sd $treat_frames_hw if $sample, a($studycontrols $covariates) vce(robust)
				
				mat b = r(table)
				local tcount: word count $treat_frames_sd $treat_frames_hw			
				forval i = 1/`tcount' {
					local b_`i' = b[1,`i'] //beta
					local se_`i' = b[2,`i'] //asymptotic s.e.
					local a_p_`i' = b[4,`i'] //asymptotic p-value
					local p_`i' = 0 //exact p-value
				}
				
				local n : di %9.0fc e(N)
				local r2_a : di %9.2f e(r2_a)
				
				*Randomization inference using RI dataset created
				preserve
    				forval batch = 1/$nbatch {
					   use if batch == `batch' using "$rianalysis/ri_data_final", clear
    				   noi di "Running ITT-frames for `v' for batch - `batch' at time: `c(current_time)'"
    				   reghdfe `v' $treat_frames_sd $treat_frames_hw if $sample, a($studycontrols $covariates) vce(robust) 
    				   mat ri_b = r(table)
    				   forval i = 1/`tcount' {
            			   local ri_b_`i' = ri_b[1,`i']
            			   *One-tailed comparison
            		       if `ri_b_`i'' >= `b_`i'' local p_`i' = `p_`i'' + 1 
    				   }  
    				} 
				restore 
				
				*Compute Exact p-values 
				forval i = 1/`tcount' {
					local p_`i' = `p_`i''/$nbatch
				}
				
				*Imputing coefficients, asymptotic and exact p-values in the output matrix
				mat `v'_f = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat `v'_f[`i',1] = round(`b_`i'',0.001)
					mat `v'_f[`i',2] = round(`se_`i'',0.001)
					mat `v'_f[`i',3] = round(`p_`i'',0.001)
				}
				
				
				*Statistical significance
				mat stars = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat stars[`i',2] = 0
					mat stars[`i',3] = 0
					if (`a_p_`i'' < 0.1)  mat stars[`i',1] = 1
					else mat stars[`i',1] = 0
					if (`a_p_`i'' < 0.05) mat stars[`i',1] = 2 
					if (`a_p_`i'' < 0.01) mat stars[`i',1] = 3 
				}
				
				*Save matrix
				matsave `v'_f, replace saving path($tables)
				*matload `v'_f, path($tables) saving
				
				*Output
				frmttable, statmat(`v'_f) store(`v'_f) substat(2) annotate(stars) asymbol(*,**,***) rtitle($rtitle_frames) ///
				addrows("Adjusted $ R^2$", `r2_a' \ N, "`n'" \ "Control Mean", `c_mean') sdec(3,3,3)
								
			********************************************************************
			***BY 2 TIMINGS
			********************************************************************
			
				*Control mean
				mean `v' if treatment_pooled==0 & $sample
				mat mean = r(table)
				local c_mean : di %9.2f mean[1,1]
				
				*Regression 
				reghdfe `v' $treat_timings_sd $treat_timings_hw if $sample, a($studycontrols $covariates) vce(robust)
				
				mat b = r(table)
				local tcount: word count $treat_timings_sd $treat_timings_hw 		
				forval i = 1/`tcount' {
					local b_`i' = b[1,`i'] //beta
					local se_`i' = b[2,`i'] //asymptotic s.e.
					local a_p_`i' = b[4,`i'] //asymptotic p-value
					local p_`i' = 0 //exact p-value
				}
				
				local n : di %9.0fc e(N)
				local r2_a : di %9.2f e(r2_a)
				
				*Randomization inference using RI dataset created
				preserve
    				forval batch = 1/$nbatch {
					   use if batch == `batch' using "$rianalysis/ri_data_final", clear
					   noi di "Running ITT-timings for `v' for batch - `batch' at time: `c(current_time)'"
    				   qui reghdfe `v' $treat_timings_sd $treat_timings_hw if $sample, a($studycontrols $covariates) vce(robust) 
    				   mat ri_b = r(table)
    				   forval i = 1/`tcount' {
            			   local ri_b_`i' = ri_b[1,`i']
            			   *One-tailed comparison
            		       if `ri_b_`i'' >= `b_`i'' local p_`i' = `p_`i'' + 1 
    				   }  
    				} 
				restore 
				
				*Compute Exact p-values 
				forval i = 1/`tcount' {
					local p_`i' = `p_`i''/$nbatch
				}
				
				*Imputing coefficients, asymptotic and exact p-values in the output matrix
				mat `v'_t = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat `v'_t[`i',1] = round(`b_`i'',0.001)
					mat `v'_t[`i',2] = round(`se_`i'',0.001)
					mat `v'_t[`i',3] = round(`p_`i'',0.001)
				}
				
				
				*Statistical significance
				mat stars = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat stars[`i',2] = 0
					mat stars[`i',3] = 0
					if (`a_p_`i'' < 0.1)  mat stars[`i',1] = 1
					else mat stars[`i',1] = 0
					if (`a_p_`i'' < 0.05) mat stars[`i',1] = 2 
					if (`a_p_`i'' < 0.01) mat stars[`i',1] = 3 
				}
				
				*Save matrix
				matsave `v'_t, replace saving path($tables)
				*matload `v'_t, path($tables) saving
				
				*Output
				frmttable, statmat(`v'_t) store(`v'_t) substat(2) annotate(stars) asymbol(*,**,***) rtitle($rtitle_timings) ///
				addrows("Adjusted $ R^2$", `r2_a' \ N, "`n'" \ "Control Mean", `c_mean') sdec(3,3,3)
								
			********************************************************************
			***BY 10 TREATMENT ARMS
			********************************************************************
			
				*Control mean
				mean `v' if treatment_pooled==0 & $sample
				mat mean = r(table)
				local c_mean : di %9.2f mean[1,1]
			
				*Regression 
				noi reghdfe `v' $treat_arms_sd $treat_arms_hw if $sample, a($studycontrols $covariates) vce(robust)
				
				mat b = r(table)
				local tcount: word count $treat_arms_sd $treat_arms_hw		
				forval i = 1/`tcount' {
					local b_`i' = b[1,`i'] //beta
					local se_`i' = b[2,`i'] //asymptotic s.e.
					local a_p_`i' = b[4,`i'] //asymptotic p-value
					local p_`i' = 0 //exact p-value
				}
				
				local n : di %9.0fc e(N)
				local r2_a : di %9.2f e(r2_a)
				
				*Randomization inference using RI dataset created
				preserve
    				forval batch = 1/$nbatch {
					   use if batch == `batch' using "$rianalysis/ri_data_final", clear
					   noi di "Running ITT-arms for `v' for batch - `batch' at time: `c(current_time)'"
    				   reghdfe `v' $treat_arms_sd $treat_arms_hw if $sample, a($studycontrols $covariates) vce(robust) 
    				   mat ri_b = r(table)
    				   forval i = 1/`tcount' {
            			   local ri_b_`i' = ri_b[1,`i']
            			   *One-tailed comparison
            		       if `ri_b_`i'' >= `b_`i'' local p_`i' = `p_`i'' + 1 
    				   }  
    				} 
				restore
				
				*Compute Exact p-values 
				forval i = 1/`tcount' {
					local p_`i' = `p_`i''/$nbatch
				}
				
				*Imputing coefficients, asymptotic and exact p-values in the output matrix
				mat `v'_a = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat `v'_a[`i',1] = round(`b_`i'',0.001)
					mat `v'_a[`i',2] = round(`se_`i'',0.001)
					mat `v'_a[`i',3] = round(`p_`i'',0.001)
				}
				
				
				*Statistical significance
				mat stars = J(`tcount',3,.)
				forval i = 1/`tcount' {
					mat stars[`i',2] = 0
					mat stars[`i',3] = 0
					if (`a_p_`i'' < 0.1)  mat stars[`i',1] = 1
					else mat stars[`i',1] = 0
					if (`a_p_`i'' < 0.05) mat stars[`i',1] = 2 
					if (`a_p_`i'' < 0.01) mat stars[`i',1] = 3 
				}
				
				*Save matrix
				matsave `v'_a, replace saving path($tables)
				*matload `v'_a, path($tables) saving

				frmttable, statmat(`v'_a) store(`v'_a) substat(2) annotate(stars) asymbol(*,**,***) rtitle($rtitle_arms) ///
				addrows("Adjusted $ R^2$", `r2_a' \ N, "`n'" \ "Control Mean", `c_mean') sdec(3,3,3)
				
		}


********************************************************************************
***COMPILE OUTPUTS
********************************************************************************

	qui foreach v of varlist $keysdhw {	
		outreg, replay(pooled_key) merge(`v'_p) store(pooled_key) 
		outreg, replay(frame_key) merge(`v'_f) store(frame_key) 
		outreg, replay(timing_key) merge(`v'_t) store(timing_key)
		outreg, replay(arm_key) merge(`v'_a) store(arm_key) 
	}

	***Globals for outputs
	global font "basefont(scriptsize) statfont(scriptsize; scriptsize; scriptsize; scriptsize) rtitlfont(scriptsize; scriptsize; scriptsize; scriptsize) ctitlfont(scriptsize) notefont(scriptsize)"
	global lines "hlines(1{0};1{0};1{0}1;1{0}1)"
		
	***Display all results
	cd "$tables"
	foreach i in pooled frame timing arm {
		outreg using itt_`i'_key, replace tex fr replay(`i'_key) $font note("") $lines ///
			ct("", "Distancing", "", "Handwashing", "" \ "", "Know", "Act", "Know", "Act") multicol(1,2,2;1,4,2)		
	}
	

log close
exit, clear
